
function [v,s,perc,sim,varargout]=see_igone( mu , sd , lb , ub , ndr )
% function [v,s,perc,sim]=see_igone( mu , sd , lb , ub , nd );
% See IG1 with mean = MU 
% and for the time being infinite variance 
% lb  (optional) lower bound of the graph (default = 0.01)  
% ub  (optional) upper bound of the graph (default = 5 )
% ndr (optional) number of draws ( default = 20000 ) 
% V=2 can change that from loc_gone.m 
% 
%
pvec=[0.005 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.995]; 
if nargin < 5 
    ndr=100000; 
    if nargin < 4 
        ub = 5 ;
        if nargin < 3 
            lb = 0.01; 
        end
    end
end
if sd > 40 
    sd=inf 
end 
[svec,vvec]=inverse_gamma_specification( mu , sd ,1 ); 
gr=linspace(lb,ub);
jj=1; 
for jj=1:length(vvec); 
    v=vvec(jj); 
    s=svec(jj); 
    pd=zeros(100,1); 
    ii=1; 
    for ii=1:100; 
        pd(ii)=pdf_igone(gr(ii),v,s) ; 
    end 
    mod=( (v/(v+1))^0.5 )*sqrt(s/v); 
    
    into.title=['IG1 Mean ',num2str(mu),' SD: ',num2str(sd), ' ; v ', num2str(v) , ' s= ',num2str(s),' mode : ',num2str(mod) ]; 
    plot(gr,pd); 
    title([' PDF ',into.title ] );
    dr=zeros(ndr,1); 
    sim=1./gamrnd(v/2,2/s,ndr,1); 
    sim=sqrt(sim);
    display( ['Mean Simulated ' num2str( mean(sim) ) '; Std simulated ' num2str( std(sim) ) ] ) ;
    pvec=floor(ndr*pvec); 
    sims=sort(sim); 

    perc=[(pvec/ndr)' sims(pvec)]; 
    nper=length(pvec); 
    cper=emptycell(nper,1); 
    
    std_sim = std( sim ); 
    
    
    ii=1; 
    for ii=1:nper; 
    temp1=sprintf('%3.3f',perc(ii,1) ); 
    temp2=sprintf('%3.3f',perc(ii,2) );
    cper(ii)={ [temp1,'  ',temp2] }; 
    end; 
    legend( cper{1}, cper{2}, cper{3},  cper{4}, cper{5}, cper{6}, cper{7}, cper{8},  cper{9} );  
    into.scale=1; 
    histplot( sim , lb:( (ub-lb)/20 ) :ub , into ); 
    
    if nargout > 4
        varargout{1}=gr';
        varargout{2}=pd;
    end
    
end;
    
